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NUMERICAL METHODS FOR THE 2-HESSIAN ELLIPTIC 
PARTIAL DIFFERENTIAL EQUATION 

BRITTANY D. FROESE, ADAM M. OBERMAN, AND TIAGO SALVADOR 


Abstract. The elliptic 2-Hessian equation is a fully nonlinear partial differen¬ 
tial equation (PDE) that is related to intrinsic curvature for three dimensional 
manifolds. We introduce two numerical methods for this PDE: the first is 
provably convergent to the viscosity solution, and the second is more accu¬ 
rate, and convergent in practice but lacks a proof. The PDE is elliptic on a 
restricted set of functions: a convexity type constraint is needed for the ellip- 
ticity of the PDE operator. Solutions with both discretizations are obtained 
using Newton’s method. Computational results are presented on a number of 
exact solutions which range in regularity from smooth to nondifferentiable and 
in shape from convex to non convex. 
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1. Introduction 


In this article we study numerical approximations of a fully nonlinear elliptic par¬ 
tial differential equation (PDE), the fc-Hessian equation. The fc-Hessian equations 
are a family of PDEs in n-dimensional space, which include the Laplace equation, 
when k = 1, and the Monge-Ampere equation, when k = n. We have already stud¬ 
ied the Dirichlet problem for the Monge-Ampere equation [FOlla, FOllb, F013]. 
Here we study the first instance of this equation which is neither the Laplacian, or 
the Monge-Ampere equation, which is the 2-Hessian equation in three dimensions, 

(1) *5*2 [w] — Uxx'Uyy T 'U'Xx'U'zz T 'Uyy’Uzz ^xy ^xz ^ yz‘ 

While the 2-Hessian equation is unfamiliar outside of Riemannian geometry and 
elliptic regularity theory, it is closely related to the scalar curvature operator, which 
provides an intrinsic curvature for a three dimensional manifold. Geometric PDEs 
have been used widely in image analysis [Sap06] . In particular, the Monge-Ampere 
equation in the context of Optimal Transportation has been used in three dimen¬ 
sional volume based image registration [HZTA04] . Scalar curvature equations have 
not yet been used in these contexts, perhaps because no effective solvers for PDEs 
involving this operator have yet been developed. The 2-Hessian operator also ap¬ 
pears in conformal mapping problems. Conformal surface mapping has been used 
for two dimensional image registration [AHTK99, GWC + 04], but does not gener¬ 
alize directly to three dimensions. Quasi-conformal maps have been used in three 
dimensions [WWJ + 07, ZG11], however these methods are still being developed. 

In this article we introduce a monotone discretization of the 2-Hessian equation 
in the three-dimensional case. A proof of convergence to the viscosity solution 
is provided. We also build a second order accurate finite difference solver which, 
while ustable if a simple iteration is used, can be modified to converge in practice. 
Numerical results are presented on solutions with varying regularity. 

We focus on the Dirichlet problem 



where D is a rectangular (three dimensional box) domain, which is natural when 
treating computationally prescribed curvature problems. (For other topologies, 
different boundary conditions need to be used. For the torus, periodic boundary 
conditions can be used. For the sphere, it is more complicated, but it is possible to 
patch together several cubic domains to obtain this topology.) 

The operator is not elliptic, unless an additional constraint is imposed, which 
corresponds loosely to the requirement that the Laplacian restricted to every two- 
dimensional plane be positive. This condition is explained in Proposition 2.6 and 
if we assume that / > 0, it reduces to 

d 2 u d 2 u 

> 0, for every orthogonal triplet of vectors ( v , w, z). 

dv dw z 
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In other words, the two dimensional Laplacian restricted to every plane is positive 
for the function u. Hence the discretizations of the operator must also enforce the 
condition above. This means that either we are working with a family of inequality 
constraints, which makes the discretization very challenging, or that we need to 
find a way to encode the constraints in the PDE. We pursue the second option for 
the monotone discretization. 

1.1. Related work on curvature equations. The 2-Hessian equation is closely 
related to a curvature PDE in three dimensions. In two dimensions there are several 
works on the evolution of curves using curvature, going back to the seminal paper 
of Osher and Sethian [OS88]. In [Obe04], a finite difference monotone scheme is 
given for the motion of level sets by mean curvature. The advantage of monotone 
discretizations is that they have a convergence proof, and convergent schemes are 
more stable and allow for faster solvers [Set95]. The surface evolver [Bra92] is a 
tool to evolve two dimensional surfaces by curvature based on the minimization of 
its energy. In [Sap06] one can find a relation between geometric PDEs and image 
analysis. For a review of the numerical methods for curvature flows see [DDE05]. 

1.2. Related work on the Monge-Ampere equation. In this paper we study 
a fully nonlinear elliptic PDE, while most of the curvature flows lead to quasi- 
linear parabolic papers. Thus, we also review some of the related work on the 
Monge-Ampere equation, a fully nonlinear elliptic PDE. For an extended review on 
numerical methods for fully nonlinear elliptic PDEs see [FGN13]. 

The Monge-Ampere equation has been exhaustively studied. Consistent schemes 
using either finite elements [Neil3, BN12] or finite differences [LR05] have been 
proposed. However, these schemes are not monotone and therefore do not fall 
within the convergence framework of Barles and Souganidis [BS91]. They require 
instead the PDE solution to be sufficiently smooth and the numerical solver to be 
well initialized. Using wide stencil discretizations, consistent monotone schemes 
were built [FOlla, FOllb], which are thus provably convergent but have limited 
accuracy due to their directional resolution. This limitation has been overcome 
recently. By introducing filtered schemes, which blend a monotone scheme with an 
accurate (but possibly unstable scheme), the authors in [F013] were able to obtain a 
provably convergent scheme with improved accuracy. Two other solutions, specific 
to particular dimensions, have been proposed as well: in the two dimensional setting 
using a mixture of finite differences and ideas from discrete geometry [BCM14] and 
in the three dimensional setting using ideas from discretizations of optimal transport 
based on power diagrams [Mir 15]. 

The Monge-Ampere problem is related to the problem of prescribed Gauss cur¬ 
vature. A numerical method for the problem of prescribed Gauss curvature can be 
found in [MO + 14]. The Gauss curvature flow is also used in image processing for 
surface fairing [EE07]. There are very few publications devoted to solving it. In the 
early work of [SG10] a quadratically constrained eigenvalue minimization problem 
is solved to obtain the solution of the 2-Hessian equation. 

1.3. Scalar curvature and the 2-Hessian equation. The Gaussian curvature 
of a two-dimensional surface is the product of the principal curvatures, k\,K 2 of 
the surface. It is an intrinsic quantity: it does not depend on the embedding of 
the surface in space. Locally, the surface can be defined as the graph of a function 
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u(x), whose gradient of the function vanishes at x. Then the Gaussian curvature 
at x is given by the determinant of the Hessian of u(x), 

det (D 2 u) = K1K2, 

which is the two dimensional Monge-Ampere operator applied to u (if the gradient 
of u does not vanish at x , additional first order terms appear). 

The sign of the Gaussian curvature characterizes the surface, and relates how 
the area of a geodesic ball in a curved Riemannian surface deviates from that of 
the standard ball in Euclidean space (larger or smaller depending on the sign). 
The uniformization theorem of complex analysis establishes the fact that every 
surface has a conformal metric of constant Gaussian curvature: the sphere, the 
Euclidean plane, or hyperbolic space. The uniformization theorem can be proved 
by several different methods. A natural method is one that solves a semi-linear 
Laplace equation for the conformal map; see [MT02, Section 8]. 

Curvature in three and higher dimensions In general dimensions, curvature is a 
tensor rather than a scalar quantity. The curvature tensor is defined by the sectional 
curvature, K ( p , a;), which is given by the Gaussian curvature of the geodesic surface 
defined by the tangent plane, p, at x. The scalar curvature (or the Ricci scalar), 
which is the trace of the curvature tensor, is the simplest curvature invariant of a 
Riemannian manifold. It can be characterized as a multiple of the average of the 
sectional curvatures. If we choose coordinates so that a three dimensional surface is 
given by the graph of a function u(x) whose gradient vanishes at x, then the scalar 
curvature is given by a constant multiple of the 2-Hessian operator: 

^ (trace(D 2 'u) 2 — trace ((D 2 u) 2 )) = K\K2 + rei«3 + K2K3 

where Ki,K2, K3 are the three principal curvatures. Again, if the gradient of u does 
not vanish at x , additional first order terms appear. However the equation above 
holds in general if we replace the principal curvatures with the eigenvalues of the 
Hessian. This leads to the 2-Hessian equation; see section 2 below. 

Since the second order terms pose the primary challenge in the solution of nonlin¬ 
ear elliptic equations, we focus on the 2-Hessian equation in this work. In a similar 
way, the Monge-Ampere equation can be related to the equation for Gauss curva¬ 
ture through the inclusion of appropriate first order terms. In [BF014] we studied 
an extension of the Monge-Ampere equation with first order nonlinear terms; in 
that case the primary challenge was the boundary conditions. 

1 . 4 . Differential geometry and /.-Hessian equations. Conformal changes of 
metric (multiplication of the metric by a positive function) have played an impor¬ 
tant role in surface theory [LP+87]. 

One of the foundational problems of Riemannian differential geometry is to gen¬ 
eralize the uniformization theorem for surfaces to higher dimensions. The general¬ 
ization of the uniformization theorem for surfaces to higher dimensional manifolds 
involves replacing constant Gauss curvature (which is a scalar in two dimensions) 
with constant scalar curvature (rather than constant tensor curvature). The result¬ 
ing problem is called, 

The Yamabe Problem Given a compact Riemannian manifold ( M , g) of di¬ 
mension n > 3, find a metric conformal to g with constant scalar curvature. 
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The solution of the Yamabe problem can be obtained by solving a nonlinear 
elliptic eigenvalue problem [Tru68]. Generalizations of the Yamabe problem to 
other curvatures result in fc-Hessian type equations [ViaOO, Via99]. 

Also of interest is 

The Calabi-Yau problem [GHJ03] Find a conformal mapping, given by U(x), 
which transforms a given metric gij to a new one given by 


9ij = exp (U)g i:j . 


The conformal mapping function U (x) satisfies a real Monge-Ampere type PDE [Yau78]. 
In certain settings (for example, the quaterionic setting), the Calabi-Yau problem 
for a manifold which is even (d = 2 n) dimensional, results in a fc-Hessian type 
equation with k = d/2 [AV10]. 

Another interesting problem where the fc-Hessian equation appears is 

The Christoffel-Minkowski Problem Find a convex hypersurface with the 
k-th symmetric function of the principal radii prescribed on its outer normals. 

It turns out that the solution of the Christoffel-Minkowski problem corresponds 
to finding convex solutions of a fc-Hessian equation on the n-sphere [GM03]. 

The 2-Hessian equation corresponds to scalar curvature, as we discuss above, 
and solving the 2-Hessian PDE (or a related one) allows for the construction of 
hyper-surfaces of prescribed curvatures, for example scalar curvature [GG02]. 

Also related are the problem of local isometric embedding of Riemannian surfaces 
in R 3 and the related Weyl problem [TW08]. 


2. Background on the equation 


In this section, we present the background analysis for the fc-Hessian equation, 
with particular focus on the 2-Hessian equation in the three dimensional case. We 
follow the review by Wang [Wan09]. 

The fc-Hessian equation can be written as 



where 1 < k < n, Sk[u] = ak(X(D 2 u)), \(D 2 u) = (Ai,..., A„) are the eigenvalues 
of the Hessian matrix D 2 u and 





is the fc-th elementary symmetric polynomial. It includes the Poisson equation 


(fc = 1) 


A u = /, 


and the Monge-Ampere equation (fc = n) 


det D 2 u = /, 


as particular cases. 

The Diric.hlet problem is given by 


(kH) 



in f2, 
on dVl. 
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Admissible functions and ellipticity. When k is even, the fc-Hessian equation lacks 
uniqueness: if u solves the k —Hessian equation, so does —u. Thus an additional 
condition is needed to ensure solution uniqueness. Moreover, when studying the 
Poisson equation it is customary to focus on the case / > 0, which is equivalent 
to look for solutions that are subharmonic since as a result a maximum principle 
holds. In the case of the Monge-Ampere equation, we impose instead the additional 
constraint that u is convex, which is required for the ellipticity of the equation. In 
either cases, it is thus necessary to restrict the solutions to an appropriate class of 
functions in order to ensure the equation has interesting properties. 

Set 

r fe = {A e M” I CTj(A) > 0,j = 1,... ,fc} . 

Tfc is a symmetric cone, meaning that any permutation of A is in T*,. When k = 1, Pi 
is the half space {A £ R" | Ai + ... + A n > 0}. When k = n, T ra is the positive cone 
T n = {A e R" | Xj > 0, j = 1,... ,n}. The result is a restriction to subharmonic 
functions for k = 1 and convex functions for k = n, as mentioned above. 

Definition 2.1. A function u £ C 2 is k—admissible if X(D 2 u) £ T^. 

Proposition 2.2. If u is k—admissible then the k—Hessian equation (kH) is (de¬ 
generate) elliptic. 

Remark 2.1. We allow the eigenvalues of u to lie in the boundary ofT k and in 
such case the k—Hessian equation may become degenerate elliptic. 

Viscosity Solutions. Well-posedness and regularity for the equation is studied in 
[CNS85]. Here we start by recalling a well posedness result. 

Definition 2.3. We say that XI C R” is (k — l)-convex if it satisfies 

<Jk- i(«) > Co > 0 on dXl 

for some positive constant cq where n = («i, ..., K n -\) denote the principal curva¬ 
tures of dXl with respect to its inner normal. 

Theorem 2.4. Assume that XI is a bounded (k— 1 )-convex domain in R™ with C 3,1 
boundary dXl, g £ C 3,1 (dXl) and f £ C 1 ’ 1 (fi) with f > fo > 0. Then there is a 
unique k-admissible solution u £ C 3,a (Xl) to the Dirichlet problem (kH) for some 
a €(0,1). 

We now recall the definition of viscosity solutions. 

Definition 2.5. A function u £ USC (ff) is a viscosity subsolution of (kH) if for 
every <fi £ C 2 (Si) D T*,, whenever, u — 4> has a local maximum at x £ fl then 

[<r k (\{D 2 (j)(x)))<f 1 if x £ Vl, 

[min (a k {\(D 2 (j>(x))) - f,u - g) < 0, if x £ dTl. 

Similarly, a function u £ LSC (fl) is a viscosity supersolution of (kH) if for every 
(f> £ C 2 (Si) nr fc , whenever, u — <j> has a local minimum at x £ It then 

ia k {X{D 2 (j)(x))) > /, if x £ XI, 

[max (a k {X(D 2 (f(x))) - f,u - g) > 0, if x £ dXI. 

Finally, we call u a viscosity solution of (kH) if u* is a viscosity subsolution and 
u * is a viscosity supersolution of (kH). 
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The equations we consider satisfy a comparison principle. 

(CP) 

Suppose (kH) has a (continuous) viscosity solution. If u £ USC (fl) is a 

subsolution and v £ LSC (fl) is a supersolution of (kH), then u < v on 17. 

The proof of this result is one of the main technical arguments in the viscosity 
solutions theory [CIL92], 

We remark that Definition 2.5 allows for discontinuous viscosity solutions. How¬ 
ever, the comparison principle (CP) does not hold in this setting. The theoretical 
details of discontinuous viscosity solutions are not well established, and are well 
beyond the scope of the present article. 

2 -Hessian equation. In this paper, we focus on the the three-dimensional case with 

k = 2 

= / 

where 

(2) S 2 M = ( 72 (A) = A 1 A 2 + A 1 A 3 + A 2 A 3 . 

The Dirichlet problem given by 

(2H) \S 2 [u] = f, inH, 

I u = g, on dfl. 

Alternative description o/T 2 . We have 

r 2 = {A £ R 3 | Ai + A 2 + A 3 > 0, ( 72 (A) > 0} . 

The following Proposition provides an alternative characterization of r 2 . 

Proposition 2 . 6 . Let 

(3) r = {A € R 3 | Ai + A 2 > 0 , Ai + A 3 > 0 , A 2 + A 3 > 0 } 

Then 

r 2 = rn{A<ER 3 |ct 2 (A)>o}. 

Proof. Proving the 3 part is straightforward. We then prove the inclusion C. 
Suppose that (Ai,A 2 ,A 3 ) £ r 2 . Without loss of generality we can assume that 
Ai < A 2 < A 3 . Thus, it is sufficient to show that Ai + A 2 > 0. Suppose that 
Ai + A 2 < 0. We consider two cases, each leading to a contradiction. 

• Ai + A 2 = 0 

We have A 1 A 2 < 0. Hence 

(7 2 (A) = A 1 A 2 + A 1 A 3 + A 2 A 3 
= A1A2 + (Ai + A 2 )A 3 
= AiA 2 
< 0 , 

contradicting the assumption ct 2 (A) > 0. 

• Ai -f A 2 < 0 
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Since Ai < A 2 , we have Ai < 0. Moreover 


02(A) > 0 Aa(Ai + A2) > — A1A2 


A3 < — 


A 1 A 2 

Ai + A2 


and 


Ai + A2 + A3 > 0 


A3 > —Ai — A2 


From the above two inequalities we get 


—Ai — A2 < — 


A 1 A 2 
Ai + A2 


which we can rewrite as 

Ai(Ai + A2) + Aj < 0 . 


Now, since Ai < 0 and Ai + A 2 < 0, the left-end side of the inequality must be 
positive and we have thus derived a contradiction. □ 


It is easy to show, using differentiation, that the function 0 2 is nondecreasing on 
the set F, which gives some insight to why the set of admissible functions is the set 
of functions where S 2 is elliptic. 

The constraint 02 (A) > 0 will be enforced automatically in our schemes by taking 
a non-negative / in the PDE (2H). Therefore it is sufficient to look at the set T 
as defined in (3). We will refer to this restriction as plane-subharmonic since it 
corresponds to u being subharmonic on every plane. 


Alternative description of the 2-Hessian operator. For a 3 x 3 matrix M, the char¬ 
acteristic polynomial is given by 

det(M) — c(M) A + trace(M)A 2 — A 3 

where c(M), the sum of the principal minors of M. is given by 

(4) c(M) = ^ (trace(M ) 2 — trace(Af 2 )) . 

If Ai, A 2 and A 3 are the eigenvalues of M then 

c(M) = A 1 A 2 + A 1 A 3 + A 2 A 3 . 

Therefore, using (2), we conclude that (1) holds, 

‘5*2['ll] — C (-D — U X x^yy 4“ ^xx^zz 4“ ^yy'^zz H xy H xz H yz • 

Linearization. The linearization of c(M) defined in (4), is given by: 

Vc(M) • N = trace(M) trace(-ZV) — trac e(MN). 

We can apply the linearization of c(M) to obtain the linearization of the 2- 
Hessian operator, S< 2 [u], for u G C 2 , 

(5) V 52 [u] • v = trac e(D 2 u) trace(D 2 i^) — trac e(D 2 uD 2 v). 

Lemma 2.7. Let 11 £ C 2 . The linearization of the 2—Hessian operator (5) is 
elliptic if u is 2-admissible. 
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Proof. Without loss of generality, we choose coordinates such that D 2 u(x) is diag¬ 
onal. We can then rewrite the linearization of the 2-Hessian operator as 

VS 2 M • v = trace (AD 2 v) 

where A = diag(A 2 + A 3 ,Ai + A 3 ,Ai + A 2 ). Hence, the linearization is elliptic if 
A is positive definite, which is true if u is 2-admissible. (It also follows directly 
from the definition of nonlinear elliptic operator (in the sense of [CIL92]) that the 
linearization is elliptic.) □ 

Remark 2.2. When the function u fails to be “strictly” 2-admissible, the lineariza¬ 
tion can be degenerate elliptic, which affects the conditioning of the linear system 
(5). When u is not 2-admissible, the linear system can be unstable. 


3. Discretization and solvers 

In this section we explain why the naive finite difference method fails in general. 
We introduce explicit, semi-implicit, and Newton solvers for the naive finite differ¬ 
ence method, which perform better by enforcing the plane-subharmonic constraint. 
This is similar to the solvers used in [BFOIO] for the Monge-Ampere equation. Then 
we introduce a discretization which is monotone and thus provably convergent. 

While the monotone discretization is less accurate, it has the advantage that it 
gives a globally consistent, monotone discretization of the operator, meaning that 
we can apply the operator to non-admissible functions. This is useful because it 
circumvents the need for special initial data, and allows for the parabolic (time- 
dependent) operator to be defined on an unconstrained class of functions. 

In addition, we could combine the monotone discretization with the naive finite 
difference discretization to obtain provably convergent, accurate filtered finite dif¬ 
ference schemes, using the ideas in [F013]. This approach combines the advantages 
of both schemes, with little additional effort. In this work, we were mainly inter¬ 
ested in comparing the performance of the two schemes, so we did not implement 
the filtered scheme. 


3.1. Naive finite difference scheme. We begin by discussing the naive finite 
difference discretization of the 2-Hessian. This is done by simply using standard 
finite differences to discretize the operator. Denote by D 2 ’ h u the discretized Hessian 
using standard finite differences on a uniform grid with grid spacing h, i.e., 


D 2 ’ h Uijk = 


PxxUijk 'P:/:yU'rfk Pxz^ijk 

't^xy'nijk 'Dyy'Uijk Pyz^ijk 

Pxz^ijk Pyz'U-'jjk Pzz'U'ijk 


where, e.g., 


'Pxx'U’ijk 
Pxy’Uijk — 


l,k ij^k T ^i,j—l,k 

’ 


Wi+l,i+l,fc + Ui- 1 J —l,fc — Wi-l,j+l,fc — Ui+l,j-l,k 

W 2 ' 

We then get the discrete version of the 2-Hessian operator S2 [u] as 
( 6 ) S£[u] = c(D 2 ’ h u) 


Since we are using centered finite differences, this discretization is consistent, and 
it is second order accurate if the solution is smooth (hence the superscript A). 
However, this scheme is not monotone due to the off-diagonal terms in the cross 







10 


BRITTANY D. FROESE, ADAM M. OBERMAN, AND TIAGO SALVADOR 


derivatives u xy , u xz and u yz . Therefore the Barles and Souganidis theory [BS91] 
does not apply and no convergence proof is available. 

3.2. Failure of the parabolic solver for the naive finite differences. In 

this section we give a simple example to illustrate that the use of the naive finite 
difference scheme (6) together with a parabolic solver fails to converge. 

The parabolic solver is given by 

(7) u n+1 =u n + dt(S 2 [u] -/). 

Consider the solution of (2H) in [0, l] 3 , given by 

w( x ) = y> /( x ) = 3. 

The iteration is initialized with the exact solution with noise from a uniform distri¬ 
bution U(— 0.01, 0.01). The result after performing two iterations with the parabolic 
solver (7) with time step dt = dx A and the initial guess are illustrated in Figure 1. 
Regardless of the time step choosen (dt = gLt 4 /10 and dt = gLt 4 /100 were also 
used), after a sufficient number of iterations the solution behaves like in the ex¬ 
ample of Figure 1, until it eventually blows up. This tells us that the instability 
of the parabolic solver is inherent from the discretization rather than being the 
result of a poorly chosen time step. This instability is due to the fact that there 
is no mechanism to pick the right solution. The discretization, being a quadratic 
equation as we will see in subsubsection 3.3.1, has two solutions: the 2-admissible 
solution we are looking for and the negative of this. 




i 


Figure 1. Failure of the parabolic solver using the naive finite 
difference scheme: section z = 0.9 of the initial guess (left) and the 
solution after 25 iterations (right). 

3.3. Solvers for the naive finite difference scheme. In this section we present 
three different sofvers for the naive finite difference scheme: a Jacobi type solver 
obtained by solving the discretization for the reference variable; a semi-implicit 
solver based on an identity that relates the Laplacian and the 2-Hessian operator; 
a Newton solver. 
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3.3.1. Jacobi solver. The accurate discretization of (2H) leads to a quadratic equa¬ 
tion for the reference variable at each grid point. To see this we introduce the 
notation 


'U j i-\-l,j,k H - 'U'i —1 ,j,k 

a 2 = 

'U , i,j-\-l,k “I - 'U'i,j—l,k 

03 = 

1 “1” ^i,j,k—l 

2 

2 

2 

'U'i-\-l,j+l,k l,j—l,k 

«5 = 

V*i—ljj-\-l,k “I - — 

a 6 = 

^i+l,j,k+l “I - l,j,k— 1 

2 

2 

2 

Ui-\,j,k+l + u i+l,j,k-l 

a s = 

u i,j + l,k+l + 

ag = 

u i,j+l,k-l + u i,j-l,k+l 

2 

2 

2 


Using (6), Sf [it] = / can be rewritten as 


^4 ( 'y ' Kl Uijk)(ai 2 Uijk) J — fijk + ^^4 y ] ( a 2p 02p+l) 


l ii<z 2 <3 


4 h 4 


P = 2 


Solving for u^k and selecting the smaller root (in order to select the locally more 
plane-subharmonic solution), we obtain 

(J) 


'Uijk — 


CL± Ci2 a 3 1 

~~3 12 


\ 


8 ^2 ( a n ~ wJ 2 + 3y^(q 2 p - «2p+i) 2 + 12 fijkh 4 . 

ii<» 2<3 p= 2 


We can now use a Jacobi iteration to find the fixed point of (J). Notice that the 
plane-subharmonic constraint is not enforced beyond the selection of the smaller 
root in (J). 


Remark 3.1. Formula (J) can be rewritten as 

Uijk = Gl + g 2 + a3 - J- \Jtrace(D 2 ’ h Uij k ) 2 + 3 (f ijk - S^’ A [u]). 

Remark 3.2. Formula (J) can also be used in a Gauss-Seidel iteration, which 
should converge faster than the Jacobi iteration. We choose not to implement it 
here since all computational results were obtained in MATLAB, which is known to 
be slow with loops. 


In order to prove the convergence of the above solver, is is sufficient to prove 
that it is monotone, which in this case is the same as showing that the value Uijk is 
a non-decreasing function of the neighboring values [Obe06]. However, this is not 
the case for (J). 


3.3.2. Semi-implicit solver. The next solver we discuss is a semi-implicit one, which 
involves solving a Laplace equation at each iteration. 

We begin with the following identity for the Laplacian in three dimensions: 

I Au| — \J (Ai/)^ — ^u 2 , x F Uyy T u 2 z T 2u X xUyy T 2 u X xUzz T 2 Uyyii zz . 

If u solves the 2-Hessian equation, then 

|Au| = \/ (Aw) 2 = y u 2 xx + u\ y + u 2 zz + 2 u 2 xy + 2 u 2 xz + 2 u 2 yz + 2/ = y/\D 2 u\ 2 + 2/. 


This leads to a semi-implicit scheme for solving the 2-Hessian equation given by 
(9) A u n+l = ^\D 2 u n \ 2 + 2f. 
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Note that if u is a 2-admissible function, then Au > 0, a condition the scheme 
enforces. 

A good initial value for the iteration is given by the solution of 

Ait 0 = yjTf. 

3.3.3. Newton solver. To solve the discretized equation 

S£[u] = f 

we can also use a damped Newton iteration 

u n+l =u n_ av n 

where 0 < a < 1. The damping parameter a is chosen at each step to ensure that 
the residual ||iS^[u n ] — /|| is decreasing. (In practice we can often take a = 1, but 
damping is sometimes needed.) The corrector v n solves the linear system 

(y u s£[u n )) v n = s£[u n ] - f. 

To setup the above equation we need the Jacobian of the scheme, which is given by 
V„S^[it] = ^2 (' D VlVl u)T> U2V2 - {T> VlV2 u)T> VlV2 

Vi ,V 2 £{x,y,z}, 

Notice that it corresponds to the discrete version of the linearization of the 2-Hessian 
equation (5). 

3.4. Monotone finite difference scheme. In this section we construct a mono¬ 
tone finite difference scheme. As we saw before, the naive approach of simply using 
standard finite differences for the terms in the Hessian matrix will not work because 
the cross derivative terms u xy , u xz and u yz are not monotone. Instead the idea is 
to use wide stencils and a rotated coordinate system in which the Hessian matrix is 
diagonal. However, this coordinate system must be found in a monotone way. This 
section is divided in four parts: first, we briefly recall why it is enough to prove that 
our scheme is consistent and degenerate elliptic (and thus monotone) to conclude 
that it is convergent; second, we extend the function oi (2) to be non-decreasing 
in R 3 ; third, we find an expression for the 2-Hesssian operator S? [it] which can 
be discretized in a monotone manner; and fourth, we present the monotone finite 
difference scheme. 

3.4.1. Convergence of consistent degenerate elliptic scheme. The convergence of our 
finite difference schemes relies, as usual, on the framework developed by Barles and 
Souganidis [BS91] and its extension in [Obe06]. 

The framework in [BS91] provides us with sufficient conditions for the conver¬ 
gence of approximation schemes to the unique viscosity solution of a PDE. 

Theorem 3.1. Consider an elliptic equation that satisfies a comparison princi¬ 
ple. A consistent, stable and monotone approximation scheme converges locally 
uniformly to the (unique) viscosity solution. 

This framework, however, does not provided a method to verify monotonicity 
and stability. The work in [Obe06] accomplishes precisely that. 

Our finite difference schemes have the form 


F[u\ = F(ui,Uj eJV (i) - Ui) 
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where N(i) is the list of neighbours of zq. We say that F is degenerate elliptic if F 
is nondecreasing in each variable. 

The following Theorem, which can be found in [Obe06], yields a simple condition 
to verify both monotonicity and stability. 

Theorem 3 . 2 . A scheme is monotone and nonexpansive in the l°° norm if and 
only if it is degenerate elliptic. 

Consequently, proving that a scheme is convergent is reduced to checking two 
conditions: consistency and degenerate ellipticity. 

3.4.2. Non-decreasing extension of the operator. In this section we find a non¬ 
decreasing extension of 02 from T to R 3 . Our ultimate goal is to build a monotone 
finite difference approximation of the 2—Hessian equation. Since we know that the 
eigenvalues of admissible solutions u belong to the set T, we are free to redefine a 2 
outside of T in order to ensure convergence. We then require an extension of a 2 
that is non-decreasing in R 3 , which is accomplished in the following Lemma. 

Lemma 3 . 3 . The function a = f o sort where sort denotes the sorting function 
and f is given by 

f(x,y,z) = x ma x(y, |x|) + x max(z, |x|) + max(j/, |x|) max(z, |x|) 
extends a 2 on T and is non-decreasing in R 3 . 

Proof. Without loss of generality, we assume that x < y < z since sorting the values 
is monotone. Moreover, we can rewrite / as 

f(x, y, z) = max (y + x, |x| + x) max (z + x, |x| + x) — x 2 . 

Suppose ( x,y,z ) € T, then we recover a 2 (x,y,z). 

Next we show that a is non-decreasing as a function of (x,y,z). We have two 
cases to consider: 

• x + y > 0 

Since x < y < z, (x,y, z) £ T and so we recover a 2 which we know to be a non¬ 
decreasing function in T. 

• x + y < 0 

Since x < y < z, x < 0. We then get a(x,y,z) = —x 2 , which is increasing since 
x < 0. 

Hence a is non-decreasing. □ 

3.4.3. Elliptic expression for the operator. In this section we build an expression 
that can be discretized in a monotone way. 

The idea is to mimic what was done for the Monge-Ampere equation in [FOllb]: 
use a matrix identity to obtain a monotone expression for the operator. 

First note that trace(M) is invariant over conjugation O t MO by orthogonal 
matrices O. Second note that trace(Af 2 ) = rn'h > m 'ii with equality when 
M is diagonal. Hence we have 


trace(M ) 2 — trac e(M 2 ) < trac e(0 T M0) 2 — y ^fO T MO) 2 t 
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and therefore 


2 c(M) 

which can be rewritten as 





(10) c(M) = min <72(diag(i?)), 

O t O=i , 

R=O t MO 

where diag(i?) = (m, r 22 , ^ 33 ) is the vector which is the diagonal of the matrix i? 
and <72 is defined by (2). Thus, we have just proved the following Lemma. 


Lemma 3.4. Let M be a 3x3 symmetric matrix and V be the set of all orthonormal 
bases 0 /K 3 : 

V = {(^,^ 2 ,^ 3 ) I e R 3 , -L Vj if i ^ j, \\vih = !} • 

Then 

(11) c(M) = min <72 1 /JM 1 / 2 , v% Mv%) . 

[VI,V2,V3)&V 


We can now use Lemma 3.4 to characterize the 2-Hessian operator of a C 2 
function by expressing it in terms of second directional derivatives of u as follows: 


( 12 ) 


= min cr 2 

(V l,V2,V3)EV 


{d 2 u d 2 u d 2 u\ 

\ dv\ ’ dv 2 ’ dv% J ' 


3.4.4. Monotone operator. We now present the monotone discretization of the 2- 
Hessian operator. 

We approximate the second derivatives using centered finite differences which 
leads to a spatial discretization with parameter h. In addition, we consider a finite 
number of possible directions v that lie on the grid, thus introducing the directional 
discretization with parameter dd. We denote the set of orthogonal basis available 
on the grid by Q. We then have 

(2 H) m S$ i [u]= min d (fD VlUl u, D V2 „ 2 u, V V3V3 u ), 

(Kl,K2,K3)£P 


where 2?,,^ is the finite difference operator for the second directional derivative in 
the direction v which lies on the finite difference grid and are given by 

V uv u(xi) = 0 9 (u(xi + his) + u(xi - his) - 2 u(xi)). 

\is\ z h z 

Depending on the direction of the vector v, this may involve a wide stencil. 

We define dO as 


wfv i 

INI 


, arccos 


w 2 v 2 

wr 


, arccos 


W 3 V 3 

M 


dO = max min max arccos 

(tUl,tU2,!«3)6V (1/1,132, 123)66 

We now define Q in more detail. Let ng denote the width of the stencil and set 
Vi = {v G Z 3 : \vi\ < 1, ||v|| ^ 0} 

and for ng >2 


We then have 


V ns = {f £ Z 3 : \vi\ < n e ,V| t |<i tv £ V ne -i} ■ 
Qn e = {(^ 1 , v 2 , v 3 ) G v 3 e : Vi _L Vj if i ± j} . 
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We will refer to the monotone schemes with respect to the number of points in 
the stencil. For instance, the monotone scheme with the stencil of length 1 (i.e., 
ng = 1) has ns + 1 = 27 points. 

Remark 3.3. Given that a is a symmetric function when implementing the mono¬ 
tone scheme we do not need to look into all the triplets in Q ne . For instance, for 
ng = 1 we only need to look for the triplets in Table 1. 


VI 

V2 

V3 

( 1 , 1 , 0 ) 

(1,-1,0) 

(0,0,1) 

( 1 , 0 , 1 ) 

(1, 0 , -1) 

(0,1,0) 

(1,0,0) 

(0,1,1) 

(0,1,-1) 

(1,0,0) 

(0,1,0) 

(0,0,1) 


Table 1. Elements of up to permutations. 


no 

i 

2 

3 

4 

5 

6 

n s 

Table 2. ns is the 

26 

nun 

98 

iber 

290 
of V c 

579 

irect 

1155 

ions a\ 

1731 

mailable in the stencil, 


i.e., ns — 



Figure 2 . Elements of Vi (blue) and elements of V 2 \ Vi (orange). 


We now give the proof of the convergence of the monotone scheme. In order to 
do that, we first need to define our scheme at the boundary. Since we choose our 
domain to be the box [0, l] 3 , the grid points are aligned with the boundary and so 
we simply have to set g at those nodes. Set 


S^[u\(x) ~ f(x), ifajefi, 
u(x) — g{x), if a; € dfl. 


(M) 


F m [u]{x ) 
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Lemma 3.5. The finite difference scheme given by (M) is degenerate elliptic. 

Proof. From the definition, the discrete second directional derivatives D vv are non¬ 
decreasing functions of the differences between neighboring values and reference 
values, Uj — it,;, where Uj is one of the neighboring values of u t in the direction 
v. The scheme (2 H) M is a nondecreasing combination of the operators min and 
a (the latter proved in Lemma 3.3 to be nondecreasing) applied to the degenerate 
elliptic terms V ul/ , and so it is also degenerate elliptic. It is also clear that u — g is 
degenerate elliptic. Hence, (M) is degenerate elliptic. □ 


Lemma 3.6. Let Xq £ LI be a reference point on the grid and (f> be a C 4 function 
that is defined in a neighborhood of the grid. Then the scheme S^A] defined in 
(2 H) m approximates (2H) with accuracy 

S^[4>\ = S 2 [<f)\ + 0(h 2 + dO). 

Proof. From a simple Taylor series computation we have 

V vv (j)(x 0 ) = ^(lEo) + 0(h 2 ). 

Using (12) we can rewrite the 2-Hessian operator as 


/ d 2 4> d 2 4> d 2 (j)\ /d 2 (f> d 2 cj) d 2 4>\ 

{v 1 ,v 2 ,v 3 )ey 2 \dv 2 ’ dn 2 ’ dv\ J 02 ’ dv 2 ’ dv 2 J 


where the Vj are orthogonal unit vectors, which may not be in the set of grid vectors 
Q. We know that by definition of d6 we have 


mm max arccos 

(i'll V 2 , "3)66 


vfis 1 

IM 


, arccos 


v% v 2 

INI 


, arccos 


V 3 V 3 

INI 


< dd. 


Let then w £ Q where the above min is attained. Then the angle between between 
each Vj and Wj is less or equal than d6 and so there is dvj such that 


Vj + dvj = 


with \\dvj\\ = O(d0). 

Now we consider the discretized problem 


<$2 [4>] = . min , cr 2 {V VlVl <. 


,V v 




A) 


Tl (fPwiWi4*:'Dw2W24*l'^'IV3W3 


<t>) 

f d 2 (f d 2 <f> d 2 f>\ 2 

a2 U w 2 'dw 2 ’dw 2 ) + °^ } 

fd 2 (f d 2 (j> d 2 (f\ 2 . 

- a 2 {dv!’dvi’dvi) +0(h +d9) 

(d 2 (f d 2 (j> d 2 (/)\ , 2 

= mm a 2 l \+0(h +d0), 

(^1,^2,^3)€V \OV 1 01*2 ®^3 J 


d 2 <f> 

dw 2 


d 2 (j) 


+ G(dO). 


where we used the fact that 
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In addition, since the set of grid vectors Q is a subset of the set of all orthogonal 
vectors V up to scaling, we find that 

min 02 ('D l y 1 ^ 1 (f> 1 D^ 2V , 2 (f>, P V3V3 (f>) > min a 2 {Pif 1 v 1 4 >,Pv 2 v 2 4 >,Pv 3 v 3 (jf) 

(Vi, V2, V3)£G (l'l,<'2,<'3)S V 


/ d 2 (j> d 2 (j) d 2 (j>\ 
(vi,iy 2 ,v 3 )ev a2 \dv\ ’ dv 2 ’ dv\ J 
Combining the two inequalities deduced above, we conclude the proof. 


0{h 2 ). 

□ 


Theorem 3.7. Suppose (2H) has a continuous viscosity solution. Letuh.de denote 
the solutions of the scheme (M) and u denote the unique viscosity solution of (2H) . 
Then, as h,dO,h/dO -A 0, Uh.de converges locally uniformly to u. 

Proof. The convergence follows from verifying consistency and degenerate ellip- 
ticity, as explained above, by the Barles and Souganidis theory [BS91]. This is 
accomplished in Lemmas 3.5 and 3.6. Notice that the PDE (2H) has a comparison 
principle (CP) as pointed out in section 2. □ 


Remark 3.4. The assumption of the existence of a continuous viscosity solution 
is required for the comparison principle. This assumption is restrictive since the 
existence result of Theorem 2.f requires smooth data, which is not the case for 
the examples considered here. In fact, continuous viscosity solutions can exist in 
much more general settings. However, a precise well-posedness result for the (weak) 
Dirichlet, problem is not presently available, and the highly technical details require 
significant additional work that is beyond the scope of the present article. 

3.5. Solvers for the monotone finite difference scheme. In this section we 
present two solvers for the monotone finite difference scheme. 

3.5.1. Parabolic solver. Using the monotone discretization S^\u\, the simplest solver 
for the 2-Hessian equation is to use the fixed point method 

(13) u n+1 =u n -a(S™[u)-f) 

which corresponds to the discrete version of the parabolic equation u t = —*$'2 [it] + / 
using a forward Euler step. The fixed point iteration will be a contraction in 
the maximum norm provided that we choose a small enough, as dictated by the 
nonlinear CFL condition [Obe06], which in this case means a = 0(h 4 ). This will 
make the solver very slow. However, since we extended oa to be degenerate elliptic 
in M 3 , this is a global solver, meaning that it will converge regardless of the initial 
guess we choose. 


3.5.2. Newton solver. As with the standard finite difference scheme, one can also 
use a (damped) Newton solver. In this case the Jacobian for the monotone dis¬ 
cretization is obtained by using Danskin’s Theorem [Bcr03] and the product rule: 


v m s 2 m m 


-2 (p v * v *u)D v , v *, 


E . 

^2eKV 2 V 3 *} ; 

^1/^2 


( 'D Vl u 1 u)'D V2V2, 5 


if Pl,* l,*u T P^^u <C 0, 
otherwise, 


where v* are the directions active in the minimum in (2 H) M , with P v * v *u < 
< T) v * v *u. Unlike the previous solver, this is a local solver, meaning that 
we need a good initial guess in order to have convergence. 
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4. Computational results 

In this section we summarize the results of a number of different examples using 
the solvers described in the previous section. These computations are performed 
on a N x TV x N grid on the cube [0, l] 3 . Unless otherwise mentioned, all solvers 
were initialized with an initial guess provided by the explicit method (J), which 
we iterate until [it”] — /1 < 10 -1 . The initial guess for the explicit method (J) 
was the exact solution with some noise from a uniform distribution. As stopping 
criteria for the Newton solver we used [S^w”] ~~ f\ < 10~ 10 where H £ {A, M}. 
Solutions were also computed using (J) and (9) with very similar results to the ones 
provided by the Newton solver being obtained. For that reason, we choose not to 
display them here. 

Remark 4.1. Notice that at points near the boundary of the domain, some values 
required by the wide stencil will not be available. For this reason and to simplify 
things, we set the exact solution at those points. However it is important to point 
out that we can use interpolation at the boundary to construct a (lower accuracy) 
stencil, thus avoiding the need to initialize with the exact solution. 

Example 4.1 (Quadratic function). We consider the case where u is a non-convex 
(but 2-admissible function) given by 

(14) u(x)= xf-^xl + 2x 2 3 , /(x) = 2. 

with x = (x\, X2, X3). In Table 3, we compare the results obtained using standard 
finite differences and the monotone schemes with different stencil sizes. For this 
example, we used the Newton solver for all schemes. 

All methods provide machine accuracy which is expected since the standard finite 
differences are exact for quadratic functions and the monotone schemes computed 
the desired directional derivative. 


Errors and order, 1 st Example 

N Standard Monotone (27-point) Monotone (99-point) Monotone (291-point) 


15 

4.441 X 10“ ib 

- 

4.441 X 10~ 1B 

- 

4.441 X 10~ lb 

- 

4.441 X 10" lb 

- 

20 

4.441 x 10 _lb 

-0.00 

8.882 x 10 -16 

-2.27 

8.882 x 10" 16 

-2.27 

6.661 x lO" 16 

-1.33 

25 

4.441 x 10 -16 

-0.00 

8.882 x 10 -16 

-0.00 

8.882 x 10 -16 

-0.00 

8.882 x 10“ 16 

-1.23 

30 

4.441 x 10 _lb 

-0.00 

1.332 x 10“ 15 

-2.14 

8.882 x 10“ 16 

-0.00 

8.882 x lO” 16 

-0.00 

35 

4.441 X 10~ 16 

-0.00 

1.332 x 10“ 15 

-0.00 

8.882 x 10" 16 

-0.00 

1.110 x 10“ 15 

-1.40 


Table 3. 

Accuracy in the l°° 

norm 

and order of 

convergence of 



the schemes for the first example using the Newton solver. 


Example 4.2 (smooth convex radial function). We consider now the case where 
u is given by 

(15) u(x) = exp 2 X °^ ^ , /(x) = (3 + 2||x-x 0 || 2 )exp(||x-Xo|| 2 ). 

The maximum errors are given in Table 4. As in the previous example we used 
the Newton solver for all schemes. 

The standard finite differences provided second order convergence, which was 
expected since the solution is smooth. The monotone schemes provided only first 
order convergence (or close to it). 
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Errors and order, 2 nd Example 


N Standard Monotone (27-point) Monotone (99-point) Monotone (291-point) 


15 

2.393 x 10“ 4 

- 

3.472 x 10“ 4 

- 

2.167 x 10“ 4 

- 

1.302 x 10~ 4 

- 

20 

1.298 x 10 -4 

2.00 

2.225 x 10 -4 

1.46 

1.518 x 10~ 4 

1.17 

1.034 x 10~ 4 

0.75 

25 

8.197 x 10“ 5 

1.97 

1.650 x 10~ 4 

1.28 

1.165 x 10~ 4 

1.13 

8.552 x 10~ 5 

0.81 

30 

5.607 x 10~ 5 

2.01 

1.346 x 10 -4 

1.08 

9.357 x 10 -5 

1.16 

7.216 x 10 -5 

0.90 

35 

4.091 x 10 -5 

1.98 

1.259 x 10“ 4 

0.42 

7.809 x 10“ 5 

1.14 

6.247 x 10~ 5 

0.91 


Table 4. Accuracy in the l°° norm and order of convergence of 
the schemes for the second example using the Newton solver. 


Example 4.3 (smooth non-convex radial function). We consider now the case 
where u is given by 
(16) 

u(x) = exp (2xf — x% + 4xg) , /(x) = 8 (l + 12x^ + + lGrrg) exp (dx\ — 2x\ + 8x%) . 

The maximum errors are given in Table 5. Once again the solutions were com¬ 
puted with a Newton solver for all schemes. 

The standard finite differences demonstrates again second order convergence. 

For the monotone schemes, the error tappers off with the grid size and we only see 
an error reduction by considering wider stencils. This tells us that the directional 
resolution error dominates the spatial resolution error. It is important to point out 
that this doesn’t contradict our theoretical results since the only thing we proved 
was that we have convergence as both h and dO go to 0, which we observe here. 


Errors and order, 3 rd Example 

N Standard Monotone (27-point) Monotone (99-point) Monotone (291-point) 


15 

3.028 x 10“ 4 

- 

3.287 x 10“ 2 

- 

1.110 x 10“ 2 

- 

5.044 x 10 -11 

- 

20 

1.669 x 10~ 4 

1.95 

3.312 x 10~ 2 

-0.02 

1.211 x 10~ 2 

-0.29 

5.617 x 10~ 3 

-0.35 

25 

1.052 x 10“ 4 

1.98 

3.305 x 10~ 2 

0.01 

1.260 x 10~ 2 

-0.17 

5.920 x 10~ 3 

-0.22 

30 

7.218 x 10 -5 

1.99 

3.311 x 10 -2 

-0.01 

1.306 x lO" 2 

-0.19 

6.396 x 10 -3 

-0.41 

35 

5.262 x 10 -5 

1.99 

3.302 x 10“ 2 

0.02 

1.339 x 10" 2 

-0.16 

6.703 x 10" 3 

-0.29 


Table 5. Accuracy in the l°° norm and order of convergence of 
the schemes for the third example using the Newton solver. 


Example 4.4 (smooth non-convex radial function). We consider another example 
of smooth radial function which is non convex but 2-admissible: 

( 17 ) u ( x ) = l°g (2 + || x || 2 )> /( x ) = ~ 4 (2+|J||2)3^ 

The maximum errors are given in Table 6. Once again the solutions were com¬ 
puted with a Newton solver, regardless of the scheme. 

As in the previous example, standard finite differences provide second order 
convergence and only with wider stencils we see a decrease in error with the grid 
size. Moreover, the monotone schemes with wider stencils also exhibit second order 
convergence (before it tappers off in the case of the 99-point stencil). 
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Errors and order, A th Example 

N Standard Monotone (27-point) Monotone (99-point) Monotone (291-point) 


15 

4.723 x 10“ 5 

- 

1.664 x 10“ 3 

- 

3.882 x 10~ 4 

- 

4.909 x 10“ 4 

- 

20 

2.564 x 10 -5 

2.00 

1.668 x 10 -3 

-0.01 

1.787 x 10~ 4 

2.54 

2.500 x 10~ 4 

2.21 

25 

1.615 x 10“ 5 

1.98 

1.674 x 10~ 3 

-0.01 

1.007 x 10~ 4 

2.46 

1.462 x 10~ 4 

2.30 

30 

1.111 x 10 -5 

1.98 

1.672 x 10 -3 

0.01 

8.617 x 10 -5 

0.82 

9.063 x 10 -5 

2.53 

35 

8.052 x 10 -6 

2.02 

1.670 x 10“ 3 

0.01 

9.620 x 10“ 5 

-0.69 

6.506 x 10" 5 

2.08 


Table 6. Accuracy in the l°° norm and order of convergence of 
the schemes for the fourth example using the Newton solver. 


Example 4.5 (non smooth convex function). We consider now the case where u 
is given by 
(18) 

«(*) = i ((||x - xo|| - 0.2)+)=, /(x) = (3 + 25||x ( xc[|p - 5||x ) Xo|| ) l(||x-„||>o. a) <*). 

The maximum errors are given in Table 7. Due to its degenerate ellipticity, the 
monotone schemes required the use of the damped Newton solver. 

Despite the lack of smoothness of the solution, the Newton solver with standard 
finite differences still converged. As for the monotone scheme, there was a signifi¬ 
cant increase in the number of iterations required: the wider the stencil, the more 
iterations required (around 10 times more iterations when compared to the Newton 
solver for the naive finite differences in the worst cases). 

For the 291-stencil, as in Example 4.3, the error tapers off, indicating that the 
directional resolution error dominates the spatial error and, again, we still see the 
convergence as both h and dd go to 0. 


Errors and order, 5 th Example 

N Standard Monotone (27-point) Monotone (99-point) Monotone (291-point) 


15 

7.580 x 10 -4 

- 

2.261 x 10~ 3 

- 

7.707 x 10~ 4 

- 

5.086 x 10- 4 

- 

20 

6.506 x 10“ 4 

0.50 

2.329 x 10~ 3 

-0.10 

7.235 x 10~ 4 

0.21 

1.924 x 10~ 4 

3.18 

25 

3.353 x 10~ 4 

2.84 

2.057 x 10 -3 

0.53 

5.871 x lO" 4 

0.89 

1.758 x 10~ 4 

0.39 

30 

3.032 x 10“ 4 

0.53 

2.156 x 10“ 3 

-0.25 

5.431 x 10~ 4 

0.41 

2.197 x 10~ 4 

-1.18 

35 

2.129 x 10 -4 

2.22 

2.018 x 10~ 3 

0.42 

5.159 x 10~ 4 

0.32 

2.351 x 10~ 4 

-0.43 


Table 7. Accuracy in the l°° norm and order of convergence of 
the schemes for the fifth example using the Newton solver. 


Example 4.6 (example with blow-up). We considered as well the case 
(19) u(x) = —^3 — ||x|j 2 , f(x) = - ~ 9+l|x| i, . 

' W V ’ (_ 3 + || x || 2 )2 

Notice that / is unbounded at the boundary point (1,1,1) and u will be singular 
at that point as well. Despite that the Newton solver still converged, but with a 
smaller rate of convergence (approximately 0.3). It is important to observe that in 
the case of the Monge-Ampere, the Newton solver failed to converge in the ana¬ 
logue example. This may be because the Monge-Ampere equation is more strongly 
nonlinear than the 2-Hessian equation. The better accuracy of the wider monotone 
schemes is explained by the fact that the exact solution is prescribed at more grid 
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points near the boundary of the (computational) domain, in particular, where u is 
singular and / is unbounded. 


Errors and order, 6 th Example 

N Standard Monotone (27-point) Monotone (99-point) Monotone (291-point) 


15 

1.104 x 10“ 3 

- 

5.627 x 10~ 3 

- 

5.600 x 10~ 4 

- 

3.026 x 10~ 4 

- 

20 

1.096 x 10“ 3 

0.02 

5.224 x 10" 3 

0.24 

4.229 x 10" 4 

0.92 

2.628 x 10~ 4 

0.46 

25 

1.054 x 10 -3 

0.17 

4.891 x 10“ 3 

0.28 

3.454 x 10“ 4 

0.87 

2.344 x 10~ 4 

0.49 

30 

1.007 x 10" 3 

0.24 

4.698 x 10~ 3 

0.21 

2.921 x 10~ 4 

0.88 

2.102 x 10~ 4 

0.58 

35 

9.621 x 10~ 4 

0.29 

4.612 x 10~ 3 

0.12 

2.538 x 10 -4 

0.89 

1.906 x 10~ 4 

0.62 


Table 8. Accuracy in the l°° norm and order of convergence of 
the schemes for the sixth example using the Newton solver. 


Example 4.7. We consider as well the example with / = 1 and g = 0 Dirichlet 
boundary conditions. No exact solution is known. In Figure 3, we illustrate some 
of the surface plots of the level sets u = c of the solution with the standard finite 
differences and monotone scheme with c £ {—0.01,-0.03,-0.07}. Note that the 
zero level set (c = 0) is the boundary of the cube [0, l ] 3 where the zero Dirichlet 
boundary conditions are prescribed. The surface plots become spheres as c de¬ 
creases, with c = —.01 being the only where there’s a tangible difference between 
the two schemes, most likely due to the expected higher accuracy from the standard 
finite differences. In Figure 4, we plot the curve u(t, t, t) with t £ [0,1] and see that 
there’s a small difference between the solutions from the standard finite differences 
and the monotone scheme. 

Example 4.8. We consider as well the example with / = 1 and g = 0 Dirichlet 
boundary conditions but with a different domain fl = fli U where 

fii = {(a;, y, 2 )el 3 :(i - 0.35 ) 2 + (y - 0.35 ) 2 + (z - 0.5 ) 2 < 0.3 2 }, 
fi 2 ={(x,y,z) SR 3 : (x - 0.65) 2 + {y - 0.65) 2 + {z - 0.5 ). 2 < 0.3 2 }. 

No exact solution is known. In Figure 5, we illustrate some of the surface plots 
of the level sets it = c of the solution with the standard finite differences and 
monotone scheme with c € {0,—0.01,—0.02,—0.03,—0.035,—0.039}. In this case 
the zero level set is not convex, with the level sets u = c becoming more convex 
with smaller values of c. In this case the difference between the standard finite 
differences and monotone scheme is even smaller than in Example 4.7, as we can 
see in Figure 6, where we plot the curve u(t,t,t ) with t £ [0,1]. 

5. Conclusions 

The 2-Hessian equation is a fully nonlinear Partial Differential Equation which 
is elliptic provided the solutions are restricted to a convex cone, which we called 
plane-subharmonic. It is natural to compare this equation with the Monge-Ampere 
PDE, which is elliptic on the cone of convex functions, and which has been studied 
numerically in previous work by two of the authors. The elliptic 2-Hessian equation 
is more challenging because the constraints for ellipticity are less restrictive. 

We gave two different discretizations for the 2-Hessian equation in the three- 
dimensional case: a naive one obtained by simply using standard finite differences 
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Figure 3. Surface plots of the level sets of the solution to Example 
4.7 on a 30 x 30 x 30 grid with the naive finite differences (left) 
and the 27-point monotone scheme (right). 



Figure 4. Plot of the curves t >->• u(t, t , t) of the solution of Ex¬ 
ample 4.7 on a 30 x 30 x 30 grid. 


to discretize the Hessian and a monotone discretization that takes advantage of 
a characterization of the operator using a matrix inequality (12). The monotone 
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Figure 5. Surface plots of the level sets of the solution to Example 
4.8 on a 30 x 30 x 30 grid with the naive finite differences (left) 
and the 27-point monotone scheme (right). 



Figure 6 . Plot of the curves t ha u(t,t,t) of the solution of Ex¬ 
ample 4.8 on a 30 x 30 x 30 grid. 

discretization is provably convergent but less accurate, because the monotone dis¬ 
cretization required the use of a wide stencil. Computational results were provided 
using exact solutions of varying regularity and shape, from smooth to non differen¬ 
tiable, and from convex to nonconvex. 
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The naive discretization failed, unless we introduced a mechanism for selecting 
the correct 2-admissible (plane-subharnronic) solution. Once this mechanism was 
introduced, experimental results on a variety of solutions demonstrated that the 
method appeared to converge. The standard finite difference discretization failed 
using a standard parabolic solver. Two alternative solvers were presented, which 
enforced the “plane-subharmonic” restriction and proved to work numerically for 
all the examples considered. Additionally, a Newton solver was also implemented, 
converging for all examples considered, even for degenerate ones or with singu¬ 
lar right-hand sides, whenever initialized with a good initial guess. For smooth 
examples, we obtained second order convergence. 

The monotone discretization, less accurate due the introduction of a directional 
resolution to make it monotone, is stable and provably convergent. Numerical ex¬ 
amples show that the directional resolution easily dominates the spacial resolution, 
a natural consequence of the three dimensional setting. 

Moreover, one could have implemented filtered schemes, previously introduced 
in [F013], which would provide schemes that are provably convergent but with 
greater accuracy than the monotone schemes. However, we did not implement 
them here, since our main goal was to compare the two different discretizations 
presented and, moreover, the accurate scheme by itself proved to be convergent for 
all the examples considered, even degenerate ones. 

The 2-Hesssian equation is related to the scalar curvature, these are equal up 
to a constant when the gradient of the function vanishes. A natural extension to 
the current work is to build schemes for the prescribed scalar curvature of a three 
dimensional graph. 

In this work, we chose the box domain since it is easier to deal with compu¬ 
tationally as the boundary conditions are easily implemented. Dealing with more 
complex boundaries requires additional work. It is challenging to obtain higher 
order at the boundary while maintaining second order directional derivatives. A 
natural approach would be a combination of filtered schemes at the boundary and 
multi-scale grids [OZ15]. Unstructured grids are another possibility, having been 
used successfully by one of the authors to solve several fully nonlinear elliptic equa¬ 
tions [Fro 15]. 
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